MODULE functions

	USE definitions
	
	IMPLICIT NONE

CONTAINS

!------------------------------------------------------------
!Routine do determine shadow zone for diagnostics
!------------------------------------------------------------
FUNCTION gradient_h(pos)

 ! 		x-1	x	x+1	
 !  					
 ! y+1	7	4	8	
 ! 					
 ! y	1	C	2	
 !					
 ! y-1	5	3	6	
 !					

  USE definitions
  
  IMPLICIT NONE
  
  INTEGER, INTENT(in), DIMENSION(2) 	:: pos
  INTEGER, DIMENSION(2)		        :: xpos, ypos, npos, ppos
  INTEGER	:: i, j, n
  REAL, DIMENSION(2,8)	:: p
  REAL :: par, gradient_h, v
  
  v = SQRT(windx*windx+windy*windy)

  ppos(1) = pos(1) + relpos(1,2)
  npos(1) = pos(1) - relpos(1,2)

  ppos(2) = pos(2) + relpos(2,2)
  npos(2) = pos(2) - relpos(2,2)

  CALL make_periodic(ppos)
  CALL make_periodic(npos)

  gradient_h = windx/v*(sandpit(ppos(1), pos(2))%slabnumber &
       - sandpit(npos(1), pos(2))%slabnumber)/slabsize &
       + windy/v*(sandpit(pos(1), ppos(2))%slabnumber &
       - sandpit(pos(1), npos(2))%slabnumber)/slabsize 

  !DO i=1,8
  !   
  !   p(1,i) = 	+(windx*relpos(1,i)+windy*relpos(2,i))/ &	
  !        (SQRT(windx*windx+windy*windy)* & 
  !        2*(REAL(relpos(1,i)*relpos(1,i)+relpos(2,i)*relpos(2,i)))&
  !        )
  !   
  !   npos = pos-relpos(:,i)
  !   mpos = pos+relpos(:,i)
  !   CALL make_periodic(npos)
  !   CALL make_periodic(mpos)
  !
  !   p(2,i) = 	p(1,i)*(-sandpit(npos(1),npos(2))%slabnumber+ &
  !        sandpit(mpos(1),mpos(2))%slabnumber)/(slabsize)
  !   IF(p(1,i)<0) p(2,i) = 0
  !   
  !ENDDO
  !
  !par = 0
  !gradient_h = 0
  !
  !DO i=1,8
  !   par = par + p(2,i)/2
  !ENDDO
  !
  !gradient_h = par
  ! WRITE(*,*) par
END FUNCTION gradient_h
!------------------------------------------------------------
!------------------------------------------------------------
!Routine do determine if a cell is in shadow zone via 2D derivate of slab hight
!Returnes atan(f'(x)) of a derivate in wind direction.
!Return value is negative of cells in shadow zone and positive for exposed cells
!------------------------------------------------------------
FUNCTION shadow_zone1(pos)

 ! 		x-1	x	x+1	
 !  					
 ! y+1	7	4	8	
 ! 					
 ! y	1	C	2	
 !					
 ! y-1	5	3	6	
 !					

	USE definitions

	IMPLICIT NONE
	
	INTEGER, INTENT(in), DIMENSION(2) 	:: pos
	INTEGER, DIMENSION(2)				:: npos
	INTEGER	:: i, j, n 
	REAL, DIMENSION(2,8)	:: p
	REAL :: par, shadow_zone1

	DO i=1,8

		p(1,i) = 	+(windx*relpos(1,i)+windy*relpos(2,i))/ &	
					(SQRT(windx*windx+windy*windy)* & 
					(REAL(relpos(1,i)*relpos(1,i)+relpos(2,i)*relpos(2,i)))&
					)

				npos = pos-relpos(:,i)

		CALL make_periodic(npos)

  p(2,i) = 	p(1,i)*p(1,i)*(+sandpit(npos(1),npos(2))%slabnumber- &
					sandpit(pos(1),pos(2))%slabnumber)/(slabsize)
		IF(p(1,i)<0) p(2,i) = 0
	ENDDO

	par = 0


	DO i=1,8
		par = par + p(2,i)/2
	ENDDO
	shadow_zone1 = 0
	IF(par>TAN(pi*15.0/180.0) .AND. par>0) THEN
		shadow_zone1 = 1
	ENDIF
		!WRITE(*,*) 180.0/pi*ATAN(par), par, shadow_zone1
END FUNCTION shadow_zone1
!------------------------------------------------------------

!------------------------------------------------------------
!Routine to determine if cell is in shadow zone.
!Calculates derivate in wind direction for every cell in neighbourhood
!and sets shadow parameter = 1 if derivate > tan(15) for one cell
!------------------------------------------------------------
FUNCTION shadow_zone2(pos)

! 		x-1	x	x+1	
 !  					
 ! y+1	7	4	8	
 ! 					
 ! y	1	C	2	
 !					
 ! y-1	5	3	6	
 !					

	USE definitions

	IMPLICIT NONE
	
	INTEGER, INTENT(in), DIMENSION(2) 	:: pos
	INTEGER, DIMENSION(2)				:: npos
	INTEGER	:: i, j 
	REAL, DIMENSION(2,8)	:: p
	REAL :: shadow_zone2

	DO i=1,8

		p(1,i) = 	-(windx*relpos(1,i)+windy*relpos(2,i))/ &	
					(SQRT(windx*windx+windy*windy)* & 
					(REAL(relpos(1,i)*relpos(1,i)+relpos(2,i)*relpos(2,i))))

				npos = pos-relpos(:,i)

		CALL make_periodic(npos)

		p(2,i) = 	p(1,i)*p(1,i)*(-REAL(sandpit(npos(1),npos(2))%slabnumber)+ &
					REAL(sandpit(pos(1),pos(2))%slabnumber))/(slabsize)
		IF(p(1,i)<0) THEN
			p(2,i) = 0
		ENDIF


	ENDDO

	shadow_zone2 = 0
	j=0
	DO i=1,8
		IF(180.0/pi*ATAN(p(2,i))>15.0) THEN
		shadow_zone2=1
		j=i
		ENDIF
	ENDDO
						
END FUNCTION shadow_zone2
!------------------------------------------------------------


SUBROUTINE make_periodic(pos)

	USE definitions
	
	IMPLICIT NONE

	INTEGER, DIMENSION(2), INTENT(inout)	:: pos

	IF(pos(1)<1) THEN
		pos(1) = pos(1) + xsize
	ENDIF
	IF(pos(1)>xsize) THEN
		pos(1) = pos(1) - xsize
	ENDIF

	IF(pos(2)<1) THEN
		pos(2) = pos(2) + ysize
	ENDIF	
	IF(pos(2)>ysize) THEN
		pos(2) = pos(2) - ysize
	ENDIF

END SUBROUTINE make_periodic

SUBROUTINE masscount(mass)

	USE definitions

	IMPLICIT NONE

	INTEGER :: i,j
	REAL, INTENT(out) ::  mass

	DO i=1,xsize
	DO j=1,ysize

	mass = mass + sandpit(i,j)%slabnumber*slabsize*slabsize

	ENDDO
	ENDDO

	mass = mass/(xsize*ysize)


END SUBROUTINE

SUBROUTINE shuffle

  USE definitions

  IMPLICIT NONE

  INTEGER :: i, j, k
  INTEGER, DIMENSION(2) :: buf
  REAL :: rand


  i = xsize*ysize
  DO k = 1,xsize*ysize
     CALL random_number(rand)
     j =INT(rand*(xsize*ysize-1))+1
     buf = position(i)%entry
     position(i)%entry = position(j)%entry
     position(j)%entry = buf
     i = i-1
  ENDDO

END SUBROUTINE shuffle
     
  !------------------------------------------------------------
  !Function to poll random cell on grid
  !------------------------------------------------------------
  SUBROUTINE poll_cell(pos, i)
    
    USE definitions
    
    IMPLICIT NONE
    
    INTEGER, INTENT(in) :: i 
    INTEGER, INTENT(out), DIMENSION(2) :: pos

    pos = position(i)%entry
    
  END SUBROUTINE poll_cell
  !------------------------------------------------------------

END MODULE functions
